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Abstract. The metastable lifetime of the square-lattice and simple-cubic-lattice 
kinetic Ising models are studied in the low-temperature limit. The simulations are 
performed using Monte Carlo with Absorbing Markov Chain algorithms to simulate 
extremely long low-temperature lifetimes. The question being addressed is at what 
temperatures the mathematically rigorous low-temperature results become valid. 
It is shown that the answer depends partly on how close the system is to fields at 
which the prefactor for the metastable decay either has a discontinuity or diverges. 



1 Introduction 

In many areas of science escape over a saddle point is an important physical 
phenomenon |E|. In some cases results can be derived in the low- noise limit 
[||. However, what is often missing is any idea of when the low-noise limit 
is applicable. This paper addresses the question of how low the noise must 
be before the metastable lifetime becomes dominated by the low-noise limit 
results. 

Here the lifetime of the metastable state of the kinetic Ising model is 
studied. The square-lattice Ising model is studied using Glauber dynamics 
[|| , as well as a dynamic derived from a coupling of a quantum spin \ system 
to a 1-dimensional phonon heat bath [^j^]. The simple-cubic dynamic Ising 
model with a Glauber dynamic is also simulated in the low-temperature limit. 

Although only specific cases are simulated, a reasonable hypothesis emerges 
as to the question of how low is a low enough temperature before the low- 
noise results are valid. The hypothesis is that it depends on how close the 
applied field is to a value where the low-temperature prefactor has a disconti- 
nuity or a divergence. As the field where the low-temperature prefactor has a 
discontinuity or divergence is approached, the low-temperature limit is seen 
at progressively lower temperatures. 
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2 Models and Methods 

The Hamiltonian of a spin-i system can be written as 7i = H. sp + 7ibath- The 
spin Hamiltonian is 

w sp = -j (!) 

nn i 

where J is the ferromagnetic (J > 0) nearest-neighbor interaction parameter 
due to the exchange coupling between spins, H is the applied external field, 
a z - is the z-component of the Pauli spin operator at site j, the first sum is 
over all nearest-neighbor (nn) pairs (4 nn for the square lattice and 6 for the 
simple-cubic lattice), and the second sum is over all N spins. 

To simulate the quantum system given by the Hamiltonian requires ex- 
plicit knowledge of the heat bath to which the spin system is coupled. In 
general, simulating the complete quantum-mechanical system to obtain the 
time dependence of the spin degrees of freedom is unnecessary. With the 
given spin Hamiltonian, 7i S p, the dynamic is determined by the generalized 
master equation |§,||J|] 

— JT [p(^) ) ^sp] m'm &m'm ^ ^ pip^nn^^mn 7m' mPi^jm' m, > 

n^m 

w m + W m > w \ - w 

7m'm = " , W m = ^ W km , (2) 

k^m 

where pit) is the time dependent density matrix of the spin system, m', n, 
k, and m denote the eigenstates of H. sp , p(t) rn i rn ={m'\p{t)\m) , and Wkm is a 
transition rate from the m-th to the fc-th eigenstate. For our spin Hamiltonian 
there are no off-diagonal components, and the generalized master equation 
becomes identical to the classical master equation of a classical spin h Ising 
system with Hamiltonian 

Hlsing = — J 2_, <JiCr j ~ H z E Gi ' 

nn i 

where Oj=±l is the classical Ising spin. The explicit dynamic for the system 
depends on the transition rates Wkm- 

Martin in 1977 ]8| used the quantum Hamiltonian of Eq. (Q). He made 
the assumptions that each spin was coupled to its own fermionic heat bath, 
and that the correlation times in the heat bath are much shorter than the 
times of interest in the spin system. He then integrated over all degrees of 
freedom of the heat bath. He found that with appropriate assumptions the 
dynamic for the classical Ising model consisted of randomly choosing a spin 
and flipping it with a probability given by 

exp(/3£W) 
exp(/3£; ol d) + exp(/3£: now ) 



Low Temperature Metastability: 



3 



where j3 = T _1 (with Boltzmann's constant set to unity), E ncw is the energy 
of the configuration with the chosen spin flipped and E \^ is the energy of 
the original spin configuration. Here qc is an attempt frequency, related to 
the microscopic coupling between the heat bath and the spin Hamiltonian. 
Note that to insure that all probabilities are between zero and one requires 
that < <7g < 1- In this paper qq — 1. This derived dynamic corresponds 
the Glauber dynamic of randomly choosing a spin, randomly choosing a 
random number r uniformly distributed between zero and one, and flipping 
the chosen spin if r < PGfiip- 

Recently with K. Park a dynamic was derived for the Ising model by 
coupling it to a phonon bath |^],^| . Here we concentrate on a one-dimensional 
phonon bath, but results for phonon baths in other dimensions have also been 
obtained ffl5|. Again the assumption is made that the correlation times in 
the heat bath are much shorter than the the times of interest in the spin 
system, and then the integration over all degrees of freedom of the heat bath 
is performed. In this case the dynamic is given by randomly choosing a spin, 
and by using a flip probability given by 



qp must be chosen so as to make all probabilities less than one. For the fields 
and temperatures simulated here, this is accomplished by choosing qp = 0.01. 

We are most interested in measuring the lifetime of a metastable state 
for the Ising model. We start with all spins up (a = +1), and apply a static 
field of strength H directed opposite to the spins (directed downward) . The 
total magnetization for the model is M = J^- Cj, so the magnetization starts 
at M = TV. We measure the time r required for the magnetization to reach 
AI = 0, since for this magnetization one has crossed the saddle point and is 
rapidly moving toward the equilibrium magnetization. The units for r used 
are Monte Carlo steps (mcs), where one mcs corresponds to one attempted 
spin flip. The physical time is proportional to the time in units of Monte 
Carlo steps per spin. However, the low-temperature limit relevant here is the 
Single-Droplet regime || (where a single nucleating droplet causes escape 
from the metastable state). Consequently, we will use mcs as our units of t. 
In the present work the measurement for r is repeated for 10 3 escapes using 
different random number sequences, and the average lifetime (t) is calculated 
from these escapes. The average lifetime in the Single-Droplet regime has the 
form 



where A is the prefactor and r(H) is the energy of the nucleating droplet. 
Note that the attempt frequency q enters this equation in a natural way, so 
that changing q will not change A. Here q will stand for qc for the Glauber 




(r) =Aexp (/JO jq 



(6) 
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dynamic or qp for the phonon dynamic. The prefactor A is a function of T 
and H, and depends on the explicit dynamic of the system. In a given field, f 
and the prefactor A at zero temperature can be obtained from the measured 
lifetimes using a linear fit to 

Tln({ T )q) =T\n(A) + R (7) 

Another way of analyzing the data pcj , if the low-temperature value of r is 
known, is to calculate an effective prefactor A e g at any finite temperature, 

A cS = (r) q exp(-/?r) . (8) 

In Eq. (^J) calculating A c g assumes that the value of F used is the zero- 
temperature limit of -T, i.e. the value of r derived in the low-noise limit. 

The dynamic of randomly choosing a spin and flipping it or not with the 
decision made by comparing a random number r with the flip probability pmp 
is the physical dynamic. Note the time is updated whether or not the spin is 
flipped. This dynamic cannot be changed without changing the physics. It can 
be implemented in a straightforward way on serial computers. However, the 
average lifetimes can be extremely long at low temperatures, making simula- 
tion in the straightforward manner unfeasible. Although the dynamic cannot 
be changed, the dynamic can be implemented on the computer in a more 
efficient fashion. One way of doing this is to use a rejection-free technique, by 
which only moves that are successful are implemented, and the time to make 
such a successful move is added to the current time. For discrete models, 
such as the Ising model, this method is called the n-fold way, and was imple- 
mented in continuous time by Bortz, Kalos, and Lebowitz in 1975 |ll[]. The 
n-fold way method can also be implemented in discrete time [Q . The n-fold 
way method can lead to exponential speed-ups in the algorithm. Unfortu- 
nately, particularly at small field values, the n-fold way exponential speedup 
is not sufficient to allow low-temperature Ising simulations to progress in a 
reasonable amount of computer time. However, the discrete-time version of 
the n-fold way can be further accelerated by realizing that the n-fold way 
algorithm uses a 1 x 1 absorbing Markov matrix to decide what will be the 
next spin configuration and to decide the time to exit from the current spin 
configuration. This Monte Carlo with absorbing Markov chain method (the 
MCAMC method) can naturally be generalized to the case of using s x s 



absorbing Markov Chains ||3|J14|Jl§ . The details of the MCAMC method as 
well as the projective dynamics method, both of which are used here, can be 
found in a recent review ||Tq| . For the simulations here, we have used s = 1, 
2, and 4 in the absorbing Markov chains in the MCAMC method. For any 
s the average lifetime, (r), and all other averages are the same as for the 
straightforward implementation of the dynamics. This is because by using 
the MCAMC algorithm the dynamic has not been changed, the dynamic has 
just been implemented on the computer in a more intelligent fashion. 
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3 Square-lattice Ising Model 

At very low temperatures (in the low-noise limit) the kinetic Ising model 
lifetimes are influenced by the discreteness of the lattice. This has allowed for 
exact calculation of the saddle point as well as the most probable route to 
the saddle point. For the square lattice with £2 — [2J/\H\\ + 1, the critical 
droplet is a square of size £2 x £2 with one row removed and a single overturned 
spin on one of the longest sides jl7|. Here |_^J denotes the integer part of x. 
Figure [j] shows one such critical droplet for £2 = 5, when | J < \H\ < i J. 
The average lifetime is then given by (r)q = A2 exp (/3-T2) with r% = 8 J^2 — 



Fig. 1. The low-temperature square-lattice nucleating droplet with £2 = 5 is shown 
2\H\(£l -£ 2 + 1). This is valid for low temperatures and for \H\ < 4J 0. 
3.1 Glauber Dynamics 

The prefactor for the square lattice with the Glauber dynamic was determined 
from absorbing Markov chain calculations to be A2 = § for £2 — 1 and A2 = | 
for £2 — 2 p4[. Recently the prefactor has been found to be A 2 — 3/[8(^2 — 1)] 
for I2 > 1 Jig. At low enough temperature these results should hold, except 
when 2J/\H is an integer. Figure |^ shows these prefactors for many values 
of H. Note that there is a discontinuity in the prefactor (but not in I^) when 
2J/\H\ is an integer. 

The exact values of r% and A2 are known, so the MCAMC data at finite 
temperatures can be compared with these predictions. Figure ^(a) shows that 
near \H\ = 2 J the predicted lines of Tln((r)) vs. T cross. If the measured 
lifetimes were to follow these expected curves, it would mean that for T > 
0.07 J the lifetime at a field of \H\ — 1.99J would be smaller than for a field 
of \H\ = 2.01J. In other words, there would be regions where the lifetime 
decreases as the field decreases. As seen in Fig. ||this does not actually occur. 
Rather, the low-temperature predictions only agree with the data at lower 
and lower temperatures as the value of \H\ approaches the value 2J/\H\ 
where the prefactor has a discontinuity. This is also demonstrated in Fig. 
|(b). Note that these MCAMC results using a Glauber dynamic are for 10 
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Fig. 2. The low-temperature square-lattice prefactors are shown vs. \H\. The solid 
lines and diamonds are the exact low-temperature results for the prefactors for 
the Glauber dynamic. The dashed curve and filled circle are the low-temperature 
prefactors for the Ising model with each spin connected to a d = 1 phonon bath. 
The curves and lines are derived analytically, while the symbols are from low- 
temperature MCAMC simulations. The prefactors are shown only for large values 
of \H\. Only three prefactors for the Glauber dynamic have been obtained at the 
discontinuity. Only prefactors for \H \ > J are shown for the phonon dynamic 



escapes from the metastable state for a 24 x 24 square lattice Ising system. 
The error estimates are from the second moment of t obtained from the 10 3 
escapes. These error estimates should be viewed as approximate, since the 
measured second moment most likely deviates substantially from the exact 
one when only 10 3 escapes are used. 

Figure || shows the results for A e g obtained at three different tempera- 
tures. Only MCAMC data points in the Single Droplet regime (as defined 
here by the ratio of the second moment of r to (r) being greater than i || ) 
are shown. It is seen that the exact low-temperature prefactor is approached 
at finite temperature more quickly for values of H that are far from the val- 
ues where A has a discontinuity. Note that in Fig.H the higher \H\ values 
start a cross-over toward the multi-droplet regime ||], and so should not be 
used in the comparison with the low-temperature prefactors. This cross-over 
depends on the system size as well as the temperature || . Figure [|(b) shows 
a portion of the data in Fig. ^(a) for £2 — 2. Again note that near a value 
of H where there is a discontinuity (here \H\ = J and 2 J) the effective pref- 
actor at finite temperature approaches the exact low-temperature prefactor 
extremely slowly. 

A recent preprint JlCj ] has performed a similar analysis of the prefactors 
for the square-lattice Ising dynamic using a generalized Becker-Doring ap- 
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Fig. 3. The low-temperature results for the square-lattice Ising model as functions 
of T. The solid lines are the exact low-temperature results, including the prefactors. 
(a) Three different fields are shown. Note that the lines cross, which would mean 
a non-monotonic behavior of the lifetime in this region, (b) Results for 1% = 3 
are shown for several values of \H\. Note that as the discontinuity at \H\=J is 
approached, the temperature required before the low-temperature results are seen 
decreases 




|H|/J |H|/J 

Fig. 4. The low-temperature results for the effective prefactor in the Single Droplet 
regime are shown as a function of \H\. The solid horizontal lines are the exact 
low-temperature prefactors. Three different temperatures are shown. Note that the 
deviation from the expected results at larger \H\ values can be mainly attributed to 
the cross-over between the Single and Multi-Droplet regimes, (a) For |J < \H\ < 
4 J. (b) Only for J < \H\ < 2 J, where l 2 =2 
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proach. The authors find qualitatively similar results to those shown in Fig. 
|]. However, unlike our results, their results do not approach the predicted 
low-temperature prefactors. They argue that this might be due to a difference 
between discrete and continuous times in the Glauber simulations. However, 
for the temperatures and field values simulated here we did not observe in 
our simulations that the prefactor value would depend on whether the time in 
the Glauber simulation is continuous or discrete. In the single-droplet regime 
relevant at low temperatures, the average lifetimes are extremely long, and 
consequently the average lifetimes should not depend in a significant way on 
whether the time in the simulation is continuous or discrete. The dynamic 
used in rcf. jl(| is not identical to the Glauber dynamic, and this difference 
between the two dynamics might account for the difference in the prefactors. 



3.2 Phonon Dynamics 

For the square-lattice Ising model coupled to a one-dimensional phonon bath 
Fig. |2| shows the exact low-temperature prefactor [Q for \H\ > J. This is 
given by the equation 

4(2|ff|-4J) + (8J-2|ff|) 
4(2|ff| -4J)(8J-2|F|) 1 ' 

for 2J < \H\ < 4J and 

|g|+2(2J-|H|) 

2 4 |ff|(2J- \H\) 1 ' 

for J < \H\ < 2 J [Q. The same nucleating droplets (except at \H\ = 2 J) are 
responsible for the decay of the metastable state for this dynamic as for the 
Glauber dynamic, so the same value of F-i is found for both dynamics (except 
at \H\ = 2J). Note that for this dynamic the prefactor is not constant when 
£2 is constant, and there is a divergence in the prefactor when \H\ = 2J 
or \H\ — 4J. The divergence in the prefactor at \H\ = 2J is due to the 
nucleating droplet at this field being two next-nearest neighbor overturned 
spins that has r — 8J rather than the value I2 = 4J for the Glauber 
dynamic. Recently, such a divergence in the prefactor has also been found 
in a model with continuous degrees of freedom ||l9f . Figure ^| shows values 
of A e g (from finite-temperature MCAMC) at two different temperatures as 
they approach the exact low-temperature results. Note that the divergence 
in the prefactor is only seen at lower temperatures. 



4 Simple-cubic lattice Ising Model 

For the kinetic Ising model with the Glauber dynamic on a simple cubic 
lattice the average lifetime in the low-temperature limit (for fixed H and 
system size) is given by (r) = A3 exp(/3r 3 )/g with 

r 3 = i2j^ _ sji 3 - 2\h\4(£ 3 - 1) + r 2 (11) 
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Fig. 5. The low-temperature results for the prefactors vs. \H\. The solid lines and 
diamond are the exact low-temperature results for the prefactors for the Glauber 
dynamic. The dashed curve and filled circle are the low-temperature prefactor for 
the Ising model with each spin connected to a d = 1 phonon bath. The curves and 
lines are derived analytically, while the symbols are from low-temperature Monte 
Carlo simulations. The effective prefactor from MCAMC data near H = 23 is 
shown for 10 3 escapes from the metastable state on a 32 x 32 lattice at T = O.fJ 
and T = 0.2J. Note that the low-temperature results are approached very slowly 
near the discontinuity 



with £ 3 = |_4J/|i2"|J +1 ifHIpOpiJ and T 2 the square-lattice nucleating droplet 
energy defined above. Equation ( |ll"l) should be valid when \H\ < 4 J and 
2J/\H\ is not an integer. The number of spins in the nucleating droplet 
is n c = £\{£ 3 — 1) + {£\ — £2 + 1)- The nucleating droplet is a cube with 
£3 overturned spins per side, with one face removed and the appropriate 
square-lattice nucleating droplet placed where the face was removed. Figure 
[| shows a nucleating droplet for J < \H\ < 2J, so £3 — 3 and £ 2 = 2, 
and the nucleating droplet has n c — 21. For this value of H the square- 
lattice nucleating droplet consists of three overturned spins in an L-shape. 
The prefactor for the simple cubic lattice with the Glauber dynamic has 
recently been found to be A 3 = [16(4 ~h + 1)(4 - for £3 > 3 @. 

The low-temperature prefactor for £3 = 2 (so 2 J < \H\ < 4 J) has not 
yet been derived, being a straightforward but lengthy calculation. The low- 
temperature prefactor for £3 — 1 (so 4J < \H\ < 6J) is easily derived to be 
^3 - 7/6. 

Figure |?](a) shows the expected low-temperature result near \H\ — 2 J 
for the simple cubic Ising model as a function of \H\. Projective dynam- 
ics [ p6|p2||23] | simulation results are also shown for T = 0.2T C . The critical 
temperature is T c ss 4.51J. If the prefactor were unity, the expected low- 
temperature result (solid line) would be independent of temperature. For 
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Fig. 6. A low-temperature nucleating droplet with £3 = 3 and £2 = 2. It consists 
of a cube of overturned spins with sides of length £3, with the spins from one face 
removed and a nucleating droplet with £2 placed on the removed face 



other prefactors, the dashed line is adjusted so the result with the prefactor 
depends on T. The simulation data are for a 32 3 Ising system using a projec- 
tive dynamic method []l6|,^2| . Note that away from a discontinuity the simula- 
tion data fall quite nicely on the expected low-temperature results including 
the predicted prefactor. However, they do not follow the expected results near 
the discontinuities (where £3 changes). Shown here is only the discontinuity 
near \H\ = 2J. This is similar to the result found in the square-lattice model. 
Namely, how low in temperature the simulation must be performed before 
the low-temperature result is valid depends strongly on how close one is to a 
discontinuity in the low-temperature prefactor. 

Note that, as seen in Eq. ( pd| ) and in Fig. ^(a), there is a discontinuity in F3 
whenever £3 changes. For all values of H where £3 changes the discontinuity is 
such that I3 decreases by 4 J. If the data were to follow these low-temperature 
results, as \H\ decreases through one of these discontinuities the average 
lifetime would decrease. This would occur not just because of a discontinuity 
in the prefactor, as in the square-lattice Ising model, but also because of a 
discontinuity in the exponential (in I3). 

Figure || shows MCAMC simulation results using the Glauber dynamic for 
the simple-cubic Ising model at \H\ = 3.5 J. For this value of H one has £3 = 2 
and £2 — 1 so the nucleating droplet is predicted to have n c — 5 overturned 
spins, and the predicted value of = 5 J. The low-temperature prefactor 
is not known for this value of H. As seen in Fig. ||, the low-temperature 
MCAMC data do not tend toward F3 = 5 J, but could reasonably tend to- 
ward 73 = 7 J. The published low-temperature predictions jl8 20 2lJ do not 
agree with our MCAMC results. The reason for this disagreement is that the 
previously reported low-temperature predictions are incomplete: they give 
accurate values for I3 from Eq. (|ll|) for some but not all values of H. This 
is because the most probable path of escape from the metastable state has 
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Fig. 7. The low-temperature results for the simple cubic lattice vs. \H\. The dashed 
line includes known exact prefactors. The symbols are projective dynamics simu- 
lation results on a 32 3 lattice at T = 0.2T C . (a) The solid lines are the exact 
low-temperature results (b) The solid lines are the exact, corrected, low- 

temperature results. The dashed line shows the known exact prefactors, which in 
this figure are only known for \H\ < | J. Note that now the exponential portion of 
the low-temperature prediction is continuous 
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Fig. 8. The low-temperature results for the simple cubic lattice as a function of T 
for \H\ = 3.5J. The lines have a slope of the best fit to the low-temperature data. 
The results at T — are the exact low-temperature results, both the uncorrected 
ones (lower point) and the corrected ones (upper point) 



12 Mark A. Novotny 



a droplet with higher energy than predicted by I3 from Eq. ( pd| ) for some 
values of H . For 3 J < \H\ < 4 J this higher energy saddle has an energy of 
28 J — 6|if I, that corresponds to three overturned spins in an L-shape, and so 
has n c = 3. For all \H\ < 4 J this droplet lies along the most probable path 
to the saddle point, but it is the highest energy point along this path only 
for 3J < \H\ < 4J, and hence is then the critical nucleating droplet. This 
scenario predicts that for |JJ| = 3.5J one should have I3 = 7 J, in agreement 
with the MCAMC results shown in Fig. |. 

Similar results hold near other values of H where £3 changes. For example, 
near \H\ = 2J the low-temperature predictions |l^ , ^o| , pll are that for n c — 5 
(2J < \H\ < 4J) one has T 3 = 40J - 10|F| and for n c = 21 (J < \H\ < 
2J) one has T 3 = 100J - 42\H\ (Fig. §). However, there is a droplet with 
15 overturned spins consisting of 2x2x3 overturned spins with a £2 = 2 
droplet of overturned spins on one of the 2x3 faces. This 15-spin droplet 
has r = 80J — 30\H\. This is the nucleating droplet for values of H when 
this r is larger than the predicted values of for n c — 5 or n c = 21. 
This occurs for | J < \H\ < 2J. At \H \ = 2 J the droplets with 5 and with 
15 overturned spins have the same energy, while at \H\ — |j the droplets 
with 15 and with 21 overturned spins have the same energy. The correct 
predicted values of ^(H) are hence continuous when this droplet with 15 
overturned spins is included. This is shown in Fig. 0(b). In general, one can 
show that all the predicted |ll|^0|,[n) discontinuities in r 3 where £3 changes 
vanish. After this study was completed another paper describing the low- 
temperature properties of the simple cubic dynamic Ising model J24| was 
brought to the author's attention. In ref. p3 it is shown that the critical 
droplet size is n c = ^3(^3 — S)(£a — 1) + (if — ^2 + 1) where (5=0 unless 
4 J + ^/16J 2 + \H\ 2 < \H\(2£ 3 - 1) in which case 5=1. For the simple cubic 
lattice, the critical value of ^(H) is thus continuous for all H, just as for the 
square lattice. 

5 Summary and Conclusions 

From long-time simulations of kinetic Ising models, the question of when the 
low-noise limit results are approached by finite-temperature simulations has 
been addressed. In particular, we have studied the limit where the system 
size and H are fixed, while the temperature is lowered. There are three main 
features that emerge: 

• The square-lattice and simple-cubic-lattice low-temperature predictions 
have the form (r) = Aexp(pr) with r(H) continuous almost everywhere, 
but with A(H) discontinuous. 

• The finite-temperature simulation results must be at a lower temperature 
to see the low-temperature predictions if the applied field is near a value 
where the low-temperature prefactor has a discontinuity or a divergence. 
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• The finite-temperature results have average lifetimes (r) which are de- 
creasing functions of \H\. This is true even if the low-temperature results 
are not always decreasing as \H\ increases. 

We have shown that these features are true for simulations in the square- 
lattice and simple-cubic-lattice kinetic Ising model with various dynamics. 
Whether or not these features of (r) can be shown to be true in general is an 
interesting question for future research. 
Acknowledgements 

Special thanks to M. Kolesik for allowing inclusion of unpublished projective 
dynamic data. Thanks to P.A. Rikvold and K. Park for many useful discus- 
sions. Partially funded by NSF DMR-0120310. Supercomputer time provided 
by the DOE through NERSC. 



References 

1. N.G. van Kampen, Stochastic Processes in Physics and Chemistry (Elsevier 
Science, Amsterdam, 1982) 

2. J. A. Bucklew, Large Deviaion Techniques in Decision, Simulation, and Estima- 
tion (Wiley, New York, 1990) 

3. R.J. Glauber: J. Math. Phys. 4, 294 (1963) 

4. K. Park, M. A. Novotny, Comp. Phys. Commun. in press; e-print 



mat/0109214 



cond 



5. K. Park, M.A. Novotny, in: Computer Simulation Studies in Condensed Matter 
Physics XIV, ed. by D.P. Landau, S.P. Lewis, H.-B. Smlttler (Springer- Verlag, 
Berlin), in press 

6. K. Blum, Density Matrix Theory and Applications, second edition (Plenum 
Press, New York, 1996), Chapter 8 

7. K. Binder, E. Stoll, Phys. Rev. Lett. 31, 47 (1973); K. Binder, H. Miiller- 
Krumbhaar, Phys. Rev. B 9, 2328 (1974) 

8. Ph. A. Martin: J. Stat. Phys. 16, 149 (1977) 

9. P.A. Rikvold, H. Tomita, S. Miyashita, S.W. Sides, Phys. Rev. E 49, 5080 
(1994); P.A. Rikvold, B.M. Gorman, in Annual Reviews of Computational 
Physics I, ed. by D. Stauffer (World Scientific, Singapor e, 1994), p. 149 



10. V.A. Shneidman, G.M. Nita, e-print |cond-mat/0201064 

11. A.B. Bortz, M.H. Kalos, J.L. Lebowitz: J. Comput. Phys. 17, 10 (1975) 

12. M.A. Novotny: Comp. in Phys. 9, 46 (1995) 

13. M.A. Novotny: Phys. Rev. Lett. 74, 1 (1995); erratum 75, 1424 (1995). 

14. M.A. Novotny: in: Computer Simulation Studies in Condensed Matter Physics 
IX, ed. by D.P. Landau, K.K. Mon, H.-B. Shiittler (Springer- Verlag, Berlin, 
Heidelberg, 1997), p. 182 



15. M. A. Novotny, Comp. Phys. Commun., in press; e-print cond-mat/010842£ 

16. M.A. Novotny: in Annual Reviews of Computational Physics IX, ed. by D. 
Stauffer (World Scientific, Singapore, 2001), p. 153 

17. E. Jordao Neves, R.H. Schon mann: Commun. M ath. Phys. 137, 209 (1991) 



18. A. Bovier, F. Manzo: e-print "ond-mat/0107376 



19. R.S. Maier, D.L. Stein, Phys. Rev. Lett. 87 022301 (2001) 



14 Mark A. Novotny 



20. D. Chen, J. Feng, M. Qian: Science in China (Series A) 40, 832 (1997) 

21. D. Chen, J. Feng, M. Qian: Science in China (Series A) 40, 1129 (1997) 

22. M. Kolesik, M.A. Novotny, P.A. Rikvold: Phys. Rev. Lett. 80, 3384 (1998) 

23. M.A. Novotny, M. Kolesik, P.A. Rikvold, Comp. Phys. Commun. 121-122, 303 
(1999) 

24. G. Ben Arous, R. Cerf, Electronic J. Prob. 1, 1 (1996) 



